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ABSTRACT 

The radiance and color of the twilight sky are calculated for single 
scattered radiation with the use of spherically symmetric models of the 
earth's atmosphere. Spherical geometry is used throughout the calculations 
with no plane parallel approximations. Refraction effects are taken into 
account through fine subdivision of the atmosphere into spherical shells 
of fixed index of refraction. Snell's law of refraction is used to calculate 
a new direction of travel each time that a photon traverses the interface 
between layers. Five different models of the atmosphere were used: a pure 

molecular scattering atmosphere; molecular atmosphere plus ozone absorption; 
three models with aerosol concentrations of 1, 3, and 10 times normal 
together with molecular scattering and ozone absorption. The results of the 
calculations are shown for various observation positions and local viewing 
angles in the solar plane for wavelengths in the range of 0.40y to 0.75y. 
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1. INTRODUCTION 

Numerous difficult problems arise in the study of twilight in the 
earth’s atmosphere whether by theoretical techniques or physical 
measurements. Each physical sunrise and sunset is unique, unrepeatable, 
and highly dynamic. The dynamic changes in intensities and spectral 
distributions (the observed color) as a function of viewing angle and 
relative position of the observer to the sun produce many problems in 
equipment construction and rapid data accumulation and analysis. Further- 
more, the observer does not know what viewing angles may yield the most 
useful information at a given time. 

The theoretical calculations for simple atmospheric models are 
repeatable, but may not reproduce all that is physically observed. This is 
to be expected. The relatively simple treatment of this article which 
assumed that the photons were scattered only once required an appreciable 
amount of computer time on the CDC 7600 high speed digital computer and 
produced a proportional amount of information to be processed. A later 
article by the authors will show how more suitable results may be obtained 
by the addition of multiple scattering effects. 

Hulburt (1953) attempted to determine quantitatively the brightness 
and color of the twilight sky. He showed the importance of ozone in 
determining the blue color of the zenith twilight sky. Rozenberg (1963) has 
examined these problems in great detail in his book. Most of his discussion 
is essentially qualitative and is based on approximate equations whose 
validity is difficult to assess. All of his results are for a single scattering 
model with no refraction. 

The most recent and comprehensive article on the subject is by Dave 
and Mateer (1968) . They show colorimetry results for five atmospheric models 
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(similar to the ones used in this article) , but do not include any 
refraction or multiple scattering effects. A similar treatment is used in 
the following calculations with the important difference that the refraction 
of the light rays is taken into account . 

2. COMPUTATIONAL TECHNIQUES 

The geometry of the problem for the solution of the observed radiance 

at a point P on the earth’s surface is illustrated in Figure 1. The 

quantity y in the figure is the minimum height of the direct solar beam above 

the earth's surface if there were no refraction (while h is the same 

m 

quantity when refraction is considered) . The radius of the earth R £ is 
taken as 6371 km in the calculation. 

The radiance I of the single scattered radiation as received by the 
observer at point P on the earth's surface is given by the equation 

W e) • + Vm (ij))f s (x) V s > (1) 

where 

0 o = solar zenith angle with respect to the vertical of the observer. 

6 = local zenith viewing angle for the observer. 

= volume scattering coefficient for Rayleigh scattering. 

P (y) = Rayleigh phase function for scattering angle y. 

6^ = volume scattering coefficient for aerosol scattering. 

P^(y) = aerosol phase function for scattering angle. y. 

F s (X) = incident solar flux upon scattering volume. This expression 

may be written as the product of *q(X)T s , where T g is fractional 
transmission for direct solar beam to the scattering point 


and Fq(X) is the extraterrestial solar flux at wavelength X. 
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Tp = fractional transmission of scattered light from the 
scattering point to the observer. 

The formulas for F g (X) and are given by 

T = exp(-/^(3 p (£) + e.U) - BU))d£) , (2) 

s o Jx j m 

and 

T = exp(-/ s (6 R (s*) + B 3 (s’) + B m (s , ))ds') (3) 

p o 

where 

= volume absorption coefficient for ozone. 

L = total distance traversed by photon from point of entrance 
into the model atmosphere to scattering point along the 
curved path. 

S = total distance from observation point to scattering point 
along the line of sight. 

The integtal in (1) was evaluated by simpson's rule from the single 
scattering contributions at the points of intersection of the solar rays with 
the line of sight from the observer. The integrals in (2) and (3) were 
reduced to sums of products of extinction coefficients and distances 
traversed through the atmospheric layers of the model. 

3. COLORIMETRY 

A system of coordinates called the tristimulus values provides a 

standard quantitative measurement of the color response of the human eye to 

a spectral distribution. The standard luminosity functions x^, y^, and 

z have been given in tabular form by many authors for intervals of 0.005 y. 

A 
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The values used in these calculations are from the tables of Wyszecki and 
Stiles (1967). 

The tristimulus values X, Y, and Z are obtained from the standard 

luminosity functions x , y , and z for a particular spectral distribution 

A A A 

of radiant flux represented by I dX from the equations 

A 

X = /i^x^dX, 

Y = /I x y A dX, (4) 

X = fI.~z.dX. 

A A 

Finally the chromaticity coordinates x, y, and z are obtained from the 
expressions 

x = X (X + Y + Z)" 1 , 

y = Y (X + Y + Z) _1 , (5) 

z = Z (X + Y + Z) _1 , 

Since the sum of x, y, and z is unity, only two of the chromaticity 
coordinates need be specified. 

The chromaticity coordinates for twilight were calculated from the 
radiance calculated from (1) and substituted in (4) and (5). The radiance 
was calculated for the wavelengths 0.40y [0.05] 0.75y. Since the radiance 
is a slowly varying function of wavelength, the calculated values were 
fitted by a cubic spline routine in order to obtain radiance values for 
any wavelength in this interval. These radiance values were weighted with 
the solar spectral irradiance distribution at the top of the atmosphere 
as given by Thekaekara (1972) and by the appropriate standard luminosity 
functions. Thus the integral in (4) was calculated and finally the 
chromaticity coordinates were obtained from (5) . 
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4. ATMOSPHERIC MODELS 

Five different atmospheric models were used for the calculations which 
are referred to as models A, B, C, D, and E. These models have the 
following characteristics: A, pure molecular scattering; B, molecular 

scattering plus ozone absorption; C, molecular scattering, ozone absorption, 
and normal aerosol concentration as tabulated by Elterman (1968) ; D, as 
model C, but with three times normal aerosol concentration; E, as model 
C, but with ten times normal aerosol concentration. 

The model atmosphere is assumed to have concetric spherical shells 

of thickness 1 km from the radius of the earth, R , to R 4- 400 km and shells 

e e 

of thickness 5 km from an altitude of 400 km to 500 km. Atmospheric 
parameters are taken as constant in each layer. The density values are 
obtained from the U. S. Standard Atmosphere 1962, Anonymous (1965). 

Atmospheric extinction coefficients for Rayleigh, ozone, and aerosol 
scattering are obtained from the tables of Elterman (1968) for altitudes up 
to 50 km. Above this altitude Rayleigh extinction coefficients are calculated 
from density values obtained from the U. S. Standard Atmosphere 1962, ozone 
extinction coefficients are assumed to be zero, and aerosol extinction 
coefficients are assumed to decrease from 50 to 80 km with the same scale 
height as from 40 to 50 km. The number density of the aerosols is held 
constant from 80 to 90 km in order to simulate an aerosol layer. Above 90 km 
the aerosol number density decreases with the same scale height already 
used. The calculated optical thickness of the entire atmosphere in the 
zenith direction as a function of wavelength is shown in Fig. 2 for each 
of the five models. 

The shadow height for a given line of sight of the observer is defined 
as the vertical height to the point along the line of sight that intersects 
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the path of a solar photon that just grazes the earth's surface. The 
shadow height is shown in Fig. 3 as a function of the zenith angle 0 of 
the line of sight taken as positive in the solar direction and negative in 
the antisolar direction. Each curve is for a particular value of 0 q which 
defines the observer's position on the earth's surface with respect to 
the solar direction (see Fig. 1). The solid curves give the shadow height 
when refraction effects are taken into account and the dashed curves are 
calculated without refraction. The dashed curves agree well with the 
results of Rosenberg (1963) . In many cases there is an important change 
in the shadow height when refraction is taken into account. These 
decreased shadow heights increase the radiance contributions from lower 
altitudes, particularly at longer wavelengths where the atmosphere is more 
transparent. This in turn changes the calculated colorimetery values. 

As another example of the importance of refraction, particularly for the 
rays that come close to the earth's surface, consider an observer at the 
terminator (0 q = 90°) viewing an unrefracted solar ray at the zenith (0 = 0°) 
that intersects the line of sight at an altitude of 2 km. The refracted 
ray intersects the line of sight instead at 0.55 km. When looking in the 
antisolar direction so that 0 = -85°, the same unrefracted and refracted 
.rays intersect the line of sight at 2.04 km and 0.51 km respectively. 

The transmission from outside the atmosphere to the perigee point is 
shown in Figs. 4 and 5 for X = 0.45p and 0.75y respectively. The transmission 
for Model A is much greater at X = 0.75p than at 0.45p because of the strong 
wavelength dependence of Rayleigh scattering. The ozone absorption is small 
at X = 0.45y and the difference between the curves for Models A and B cannot 
be shown on the scale of the figure. The transmission is always less at a 
shorter wavelength than at a longer when a given model and perigee height 
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are compared. The greatest changes between these figures occur in those 
regions of Fig. 5 with relatively high transmission and a perigee height 
near the earth's surface. 

5 . CALCULATED RADIANCE 

The single scattered radiance was calculated for the following observer 
positions: 0 q = 80°, 85°, 90° [2°] 108°. At each observer position 

calculations were made for the following zenith viewing angles in the solar 
plane: 0 = 90° , 85°, 80°, 75°, 70° [10° ] -70°, -75° , -80°, -85°, where 

positive values of 0 refer to viewing angles in the solar direction from 
the zenith while negative values are in the antisolar direction. All 
calculations we re performed for 8 wavelengths: 0.40y [0.05] 0.75y. 

The radiance calculated for Model A is shown in Figs. 6 and 7 as a 
function of wavelength for 0 q = 90° and various viewing angles 0. All 
calculated radiance values assume an incident solar flux of unity perpendicular 
to the incoming solar beam at each wavelength. Later these values are 
weighted by the solar distribution for the calculation of the colorimetry 
results. In the solar direction. Figure 6 shows that the greatest change 
with wavelength occurs when 8 = 85°, the large decrease in the blue caused 
by the increased Rayleigh scattering at shorter wavelengths. As the viewing 
angle increases from 85° to the zenith, the radiance in the blue increases 
at first and then decreases. There are three competing effects as 0 decreases: 
first, an increase due to less absorption of solar radiation along the ray to 
the scattering center; second, the number of scattering centers along the 
line of sight decreases by the secant effect; third, the absorption decreases 
from the scattering center to the observer. The radiance in the antisolar 
directions is shown in Fig. 7. 



- 8 - 


The radiance in the solar directions for Model B is given in Fig. 8 
for 0 q = 90°. The radiance for a fixed value of 0 now shows a pronounced 
minimum near X = 0.6p due to the ozone absorption. 

The radiance for Model C when 0^ = 90° is given in Fig. 9 for viewing 
angles in the solar direction. Rayleigh scattering, ozone absorption, 
and aerosol scattering are all included in this model. The dip in the 
radiance curves in the middle of the visible spectrum due to ozone 
absorption is clearly visible here. In most cases the radiance is less 
for Model C than for Model A except at long wavelengths and for viewing 
angles near the horizon. In this case the radiance is larger due to the 
strong forward scattering of the sunlight by the aerosols and the additional 
scattering centers along the line of sight. 

A more complete understanding of these curves can be obtained by a 
study of the differential radiance curves shown in Fig. 10. The differ- 
ential radiance per kilometer (measured along the slant path) is given as 
a function of the vertical height to the scattering volume measured along 
the earth's local radius through the infinitesimal volume. Each curve for 
a particular value of 0 increases to a maximum value; beyond the maximum 
it decreases approximately as an exponential. The initial increase in the 
radiance occurs as the transmission of the solar radiation to the 
scattering point increases with height. The differential radiance 
eventually decreases as the number of scattering centers of all types 
decreases with height. At altitudes above 60 km these curves are identical 
on the scale of this figure for the same zenith angle in the solar and 
antisolar directions (e.g. 0 = 70° and 0 = -70°) . The curves for an 
observer at the terminator have this symmetry at altitudes sufficiently 
high that there is virtually no loss in the solar radiation to the scattering 
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point on either line of sight. At lower altitudes the differential 
radiance is lower in the antisolar direction because of the additional 
absorption of the solar beam before it reaches the scattering point 
(e.g., the curve for 0 = -85° is considerably lower below an altitude of 
40 km than that for 0 = 85°) . Since the radiance observed at a given 
angle of view is proportional to the area under the corresponding 
differential radiance curve in Fig. 10, the variation in the radiance 
curves with 0 can be understood in terms of' the changes in the absorption 
of the solar beam up to the point of scattering, the dependence of the 
scattering on the phase function of the aerosols and the Rayleigh 
scattering centers, the decrease of the scattering centers with height, 
and the secant variation of the number of scattering centers along the 
line of sight. 

The radiance for Model B and 0 = 96° is shown in Fig. 11. The 

o 

radiance for Model C and 0^ = 96° is given in Figs. 12 and 13 for the 
solar and antisolar directions respectively. In general the difference 
in radiance between Models B and C is greater at long wavelengths than 
at short. This is largely due to the much greater difference between 
Models B and C in the transmission through the atmosphere up to the 
scattering point at longer wavelengths (Fig. 5) than at shorter (Fig. 4). 
The radiance in the antisolar direction shows little dependence on 0 
until 0 = -40°. The radiance rapidly decreases to very small values for 
larger zenith angles. 

The differential radiance for 0 q = 96° is given in Fig. 14. The 
differential radiance at the zenith reaches a maximum at about 62 km. 

The maxima for 0 = 40° and 0 = -40° occur approximately at 57 and 67 km 
respectively. The curves for 0=0°, 40°, and -40° are nearly the same on 
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th e scale of the figure above 75 km. Below 90 km the differential radiance 
for 0 = -70° is very much less than that for 8 = 70° due to the much 
greater atmospheric path that the solar rays must traverse to reach the 
scattering point in the former case. The maxima for these two angles are 
at about 50 and 82 km respectively. The maximum for 6 = 85° is at about 
35 km. No single scattered radiation can be observed at 0 = -85° up to 
the altitudes considered here. 

6. CHROMATICITY 

The chromaticity was calculated for each model and solar zenith angle 

by the method described in Section 3. The results are shown in a series of 

chromaticity diagrams in which the x and y chromaticity coordinates are 

plotted for each viewing angle in the solar and antisolar half planes. For 

further reference the chromaticity coordinates are given in Tables 1 to 4 

for each of the five models, for various values of 0, and for 0 = 90°, 

o 

92°, 94°, and 96°. 

Model A . Molecular atmosphere without ozone . The chromaticity diagrams 
for Model A are shown in Figs. 15 and 16 for 0 q = 90°, 92°, and 96°. The 
regions traditionally assigned to the different colors are indicated in 
Fig. 15 for guidance. The spectral purity isopleths are indicated in all 
the figures . When the spectral purity is less than about 5% the color is 
usually assumed to be a white or gray without discernable color. The 
numbers beside the curve indicate the zenith angle of observation (positive 
for the solar half plane and negative for the antisolar half plane) . The 
colors for 0 q = 90° are all yellow or an essentially white color near the 
zenith. At 0 Q = 92° and 96° the colors near the horizon have slightly greater 
spectral purity. A portion of these latter curves is also in the orange 
portion of the spectrum. The zenith sky is not blue in this model. The 
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actual colors of the twilight sky are not well represented by this model, 
except for the predominantly yellow-orange color just above the solar 
horizon. 

For comparison, results were obtained in some cases both with and 
without consideration of the effects of refraction in the atmosphere. For 
example, at 0 q = 96° and 6=0°, x = 0.351 and y = 0.346 when refractive 
effects are taken into account, but x = 0.316 and y = 0.327 without 
refraction. In all cases calculated the x and y values were less without 
refraction than the same values obtained with refraction. Thus values 
calculated without refraction tend to be too far down on the chromaticity 
diagram and in many cases would show blue or purple colors that do not 
exist when refraction effects are included. This may explain some of the 
differences between our results and those of Dave and Mateer (1968) . Other 
differences include the use of somewhat different atmospheric models and 
the extension of the atmosphere to much higher altitudes (500 km) in the 
present calculations. 

Model B . Molecular atmosphere with ozone . The results when 0.35 cm of 
ozone are added to the molecular atmosphere are shown in Figs. 17 and 18. 

The main change brought about by the addition of the ozone to the model is 
the development of blue color over a large area of the sky surrounding 
the zenith. This results from the ozone absorption by the Chappius bands 
in the vicinity of A = 0.6y. The spectral purity of the zenith increases 
as 0 q increases from 92° to 96° . A portion of the sky is purple of a low 
spectral purity in the antisolar half plane in some of these cases . The 
chromaticity coordinates vary hardly at all near the zenith as 0 changes. 

The actual turning point of the chromaticity curve occurs in the range 
from 0 = 20 to 30° for 0 q = 94° and 96°. The turning point on the figures 



- 12 - 


has usually been marked 0°, since the chromaticity is the same at 9 = 0° 
and 30° on these figures. 

The blue color of the zenith sky as well as a region of purple sky 
in the antisolar half plane are features which are obtained from this model 
which can not be explained by a purely Rayleigh scattering atmosphere. 

Model C . Molecular atmosphere , ozone , normal aerosol amount . The 

results for this model of the normal atmosphere are shown in Figs. 19 and 20. 

The blue region of the sky around the zenith has a higher spectral purity 

when aerosols are added to the model. Dave and Mateer (1968) have pointed 

out that a normal aerosol amount contributes to an increased blueness of the 

zenith sky. The spectral purity is somewhat higher in our calculations 

than in theirs and becomes almost 40% near the zenith when 0 = 96°. The 

o 

greater increase in absorption at blue than at red wavelengths for the solar 
radiation up to the scattering point when Models B and C are compared (Figs. 

4 and 5) is the reason for the deeper blue color at the zenith when aerosols 
are present. Some purple light of low spectral purity is predicted for 
some angles of observation. 

The color of the sky around the solar horizon is correctly predicted 
as an orange-red of high spectral purity. The single scattering model 
predicts orange colors near the antisolar horizon as 6 q increases in 
disagreement with observation. Obviously multiple scattering effects are 
most important in this region of the sky. The chromaticity curve for 
0 q = 94° is not shown here since it is very nearly the same as for 0 Q = 96 c . 
The main difference is that the point 0 = 80° corresponds to a yellow-orange 
of spectral purity 18% when 0 q = 94°, but corresponds to a nearly white 
color when 0 = 96° . 
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Model D . Molecular atmosphere , ozone , three times normal aerosol 
amount . The chromaticity diagram when the atmosphere contains three times 
the normal aerosol amount at all altitudes is shown in Fig. 21. At 
0 q = 92° the portion of the sky near the solar horizon has a higher spectral 
purity and tends to appear a deeper red with the greater aerosol amount. 

The additional aerosols also make the sky around the zenith a blue of 
greater spectral purity and this color extends over a larger region of 
the sky. The more vivid colors which are observed at twilight during 
periods of high aerosol concentration in the atmosphere are confirmed by 
this model. 

Model E . Molecular atmosphere , ozone , ten times normal aerosol 
amount . The chromaticity diagram when the atmosphere contains ten times 
the normal aerosol amount at all altitudes is shown in Fig. 22. This 
model represents an extreme condition with the optical thickness of the 
atmosphere in a vertical direction between 2 and 3.5 (Fig. 2) . The region 
of the sky around the solar horizon with red and orange colors of high 
spectral purity continues to expand as the aerosol concentration increases 
up to this amount. For example, at 0 q = 92° the spectral purity is greater 
than 20% from the solar horizon to 0 - 55° for Model E; the corresponding 
values of 0 for Models C and D are 78° and 76° . Furthermore the deep blue 
sky near the zenith has disappeared in Model E which has only a gray blue 
region within 30° of the zenith. These same features are also evident at 
0 Q = 96°, although there is a larger region of blue sky around the zenith 
than at 0 Q = 92°. For large aerosol amounts the model indicates vivid 
colors near the solar horizon, but colors of reduced spectral purity near 
the zenith. Multiple scattering effects would be especially important 
for the larger optical depths of this model. 
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7. CONCLUSIONS 

Molecular scattering alone cannot explain the colors of the twilight 
sky. The addition of ozone absorption is essential to a model in order to 
reproduce the blue of the zenith sky. Some purple colors of low spectral 
purity are also obtained from this model. The deep blue of the zenith sky 
and the red and orange colors of high spectral purity in the region of the 
sky around the solar horizon are obtained only when a realistic model is 
used which combines molecular, ozone, and aerosol scattering and 
absorption. When the aerosol amount is three times normal the blue of the 
zenith sky becomes deeper and extends over a greater area. However, when 
the aerosol amount is ten times normal the spectral purity of the blue of 
the zenith sky decreases appreciably. Larger regions of red and yellow colors 
of higher spectral purity around the solar horizon are predicted when the 
aerosol amount is greater than normal. 

Refraction effects have been included in the present calculations. 

The greatest limitation of the present results is that they include only 
single scattering effects. We plan to use a Monte Carlo technique to extend 
these results to include all orders of multiple scattering. 
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Legends for Figures 

Fig. 1. Geometry of direct solar flux to scattering point on line of sight 
and then to observer on the surface of the earth. 

Fig. 2. Optical thickness through the entire atmosphere in a vertical direction 
as a function of wavelength for Models A (molecular absorption) , B (ozone 
added) , C (normal aerosol amount added) , D (three times normal aerosol 
amount added) , and E (ten times normal aerosol amount added) . 

Fig. 3. Shadow height in km as a function of 8, the zenith angle of observation. 

The solar horizon is at the right and the antisolar horizon at the left. 

The angles marked on the curves are 9 as defined in Fig. 1. The dashed curves 

o 

are calculated without refraction, while the solid curves include the effect 
of the refraction of the light rays by the earth's atmosphere. 

Fig. 4. The transmission of the solar radiation to the zenith direction as 
seen at the terminator as a function of perigee height of the solar ray. 

Curves are shown for Models A, C, D, and E (the curve for Model B is the same 

as that for Model A on the scale of the figure) and for X = 0.45y. 

/ 

Fig. 5. Same as Fig. 4 except X = 0.75y. 

Fig. 6. Radiance as a function of wavelength (microns) for Model A and 0 q = 90°. 
The zenith angle of the direction of observation in the solar half plane is 
marked on each curve. 

Fig. 7. Radiance as a function of wavelength (microns) for Model A and 0^ = 90°. 
The zenith angle of the direction of observation in the antisolar half plane is 
marked on each curve. 

Fig. 8. Radiance for Model B and 0 q = 90°. 

Fig. 9. Radiance for Model C and 0 q = 90°. 



- d <L - 


Fig. 10. The differential radiance per kilometer (measured along the slant 
path) as a function of the vertical height to the scattering volume measured 
along the earth’s local radius through the infinitesimal volume. The results 
are for Model C, 0^ = 90°, and X = 0.4y. Results are shown for 0 = 85°, 

70°, 40°, 0°, -40°, -70°, -85° (negative values are in the antisolar half 
plane) . 

Fig. 11. Radiance for Model B and 0 = 96°. 

Fig. 12. Radiance for Model C and 0^ = 96° for solar half plane. 

Fig. 13. Radiance for Model C and 0^ = 96° for antisolar half plane. 

Fig. 14. Differential radiance for Model C, 0^ = 96°, and X = 0.4y. 

See caption to Fig. 10. 

Fig. 15. Chromaticity diagram for Model A and 0 q = 90°. The numbers on 
the curve indicate the value of 0 (positive in the solar half plane and 
negative in the antisolar half plane) . The spectral purity in percent is 
marked on the isopleths. The coordinates x and y are the chromaticity 
coordinates (see (5)). The generally accepted regions for the various colors 
are indicated on the diagram. 

Fig. 16. Chromaticity diagram for Model A and 0 q = 92° and 96°. 

Fig. 17. Chromaticity diagram for Model B and 0 q = 90° and 92°. 

Fig. 18. Chromaticity diagram for Model B and 0 q = 94° and 96°. 

Fig. 19. Chromaticity diagram for Model C and 0^ = 90° and 92°. 

Fig. 20. Chromaticity diagram for Model C and 0 q = 96°. 

Fig. 21. Chromaticity diagram for Model D and 0 q = 92° and 96°. 

Fig. 22. Chromaticity diagram for Model E and 0 q = 92° and 96°. 
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